% This code produces results in Table 5, 6, and 7. It also produces Figure
% 5.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
clc

addpath(genpath('LogLinearization'))
addpath(genpath('Parameter'))
addpath(genpath('Simulation'))
addpath(genpath('SubRiskPremia'))
addpath(genpath('Utility'))

addpath(genpath('jplv7'))

%% ** 0.Parameters

%% 0.1 Fundamental Parameters
nu_q = 0.018;
l_q = 372;%1927Q1--2019Q4
l_q_burn = 220;%1872Q1--1926Q4

l_a = l_q/4;

y_end = 2019;

%% 0.2 Moment Parameters

Mom_RowNames = {'$\E(\Delta c)$','$\sigma(\Delta c)$','$\E(\Delta d)$','$\sigma(\Delta d)$', ...
                '$\sigma(\tilde{\mu})$', '$\rho(\tilde{\mu})$',...
                '$\sigma(\tilde{\mu_d})$', '$\rho(\tilde{\mu_d})$',...
                '$\sigma(\tilde{\mu_r})$', '$\rho(\tilde{\mu_r})$',...
                '$corr(\tilde{\mu} , \tilde{\mu_d})$', '$corr(\tilde{\mu},\tilde{\mu_r})$', '$corr(\tilde{\mu_d},\tilde{\mu_r})$', ...
                '$\E(R_m - R_f)$', '$\sigma(R_m - R_f)$' ,'$SR(R_m - R_f)$' , ...
                '$\E(p - d)$', '$\sigma(p - d)$', '$\rho(p - d)$', ...
                '$\E(r_f)$', '$\sigma(r_f)$'};
            
MomVar_Names = {'cg_yoy_mean' , 'cg_yoy_std' , 'dg_yoy_mean' , 'dg_yoy_std' ,...
                'tmu_std' , 'tmu_rho' , 'tmud_std' , 'tmud_rho' , 'tmur_std' , 'tmur_rho' ,...
                'tmu_tmud_corr' , 'tmu_tmur_corr' , 'tmud_tmur_corr' ,...
                'EXRET_mean' , 'EXRET_std' , 'EXRET_sr' ,...
                'pd_mean' , 'pd_std' , 'pd_rho' ,...
                'rf_mean' , 'rf_std'};

type_avg = 1;% 1 for median, 2 for mean


%% ** 1. Data
Data_sum_file = 'Data/Data_Sum.xlsx';
Data_sum_q = readtable(Data_sum_file,'ReadVariableNames',1,'sheet','Quarterly','basic',true);
Data_sum_a = readtable(Data_sum_file,'ReadVariableNames',1,'sheet','Annual','basic',true);
Data_sum_pred = readtable(Data_sum_file,'ReadVariableNames',1,'sheet','Pred','basic',true);

%% 1.1 Calculate data moments
%-------------------
% Endowment
%-------------------
DATA.cg_yoy = Data_sum_q.cg(end - l_q + 1 : end);%Note that we calculate cg as 4-quarter change
DATA.dg_yoy = Data_sum_q.dg(end - l_q + 1 : end);%Note that we calculate dg as 4-quarter change

%-------------------
% Price
%-------------------
DATA.rf = Data_sum_q.tb3mreal_ea(end - l_q + 1 : end);
DATA.EXRET = Data_sum_q.Ret(end - l_q + 1 : end) - (1 + Data_sum_q.TB3m_al(end - l_q : end - 1)).^(1/4) + 1;
DATA.pd_ttm = Data_sum_pred.pd_2(end - l_q + 1 : end);

%-------------------
% State variables
%-------------------
DATA.tmu = Data_sum_pred.expc(end - l_q + 1 : end);
DATA.tmud = Data_sum_pred.expd(end - l_q + 1 : end);
DATA.tmur = Data_sum_pred.expr(end - l_q + 1 : end);


DataVar_Names = {'cg_yoy' , 'dg_yoy', 'rf' , 'EXRET' , 'pd_ttm' , 'tmu' , 'tmud' , 'tmur'};
Mom_DATA = Mom_Annual(struct2cell(DATA) , type_avg , DataVar_Names , MomVar_Names , Mom_RowNames , 'Data');

%% 1.2 Estimate alpha (Internet appendix Table B.1)
% Estimate alpha
alpha_matrix = nan(3 , 2);

dc_q = Data_sum_q.d(end - l_q + 1 : end) - Data_sum_q.c(end - l_q + 1 : end);

% Full sample estimates of alpha with different bias adjustment
dc_a = Data_sum_a.d(end - l_a + 1 : end) - Data_sum_a.c(end - l_a + 1 : end);
[dc_a_ar , dc_a_ar_adj] = VAR1_bsadj(dc_a , 10000 , 1 , 0);
alpha_matrix(1 , 1) = dc_a_ar(2)^0.25;
alpha_matrix(2 , 1) = dc_a_ar_adj(2)^0.25;
alpha_matrix(3 , 1) = AR1_Adj(dc_a_ar(2) , l_a , 2)^0.25;

% Use the sample period from Bansal et al. (2007)
dc_a_idx = Data_sum_a.d(end - (y_end - 1929) : end - (y_end - 2001)) - Data_sum_a.c(end - (y_end - 1929) : end - (y_end - 2001));%1929--2001
[dc_ar_idx , dc_ar_idx_adj] = VAR1_bsadj(dc_a_idx , 10000 , 1 , 0);
alpha_matrix(1 , 2) = dc_ar_idx(2)^0.25;
alpha_matrix(2 , 2) = dc_ar_idx_adj(2)^0.25;
alpha_matrix(3 , 2) = AR1_Adj(dc_ar_idx(2) , length(dc_a_idx) , 2)^0.25;

% Other parameters
disp('Mean and S.D. of Quarterly d-c')
disp(mean(dc_q))
disp(std(dc_q))
disp('Mean and S.D. of Annual d-c')
disp(mean(dc_a))
disp(std(dc_a))


%% ** 2. Model Simulations

%% 2.1 Parameters
c_0 = Data_sum_q.c(nonnan_first(Data_sum_q.c));
N_sim = 50000;
phi = 0.99;
psi = 1.5;
node_n = [3 , 3];%Pre-calculated with a seperate function

SimVar_Names_EIS_1 = {'cg_yoy' , 'dg_yoy', 'rf' , 'Rf' , 'Ret' , 'ret' , 'ret_MS' , 'EXRET' , 'exret' , 'pd_ttm' , 'tmu' , 'tmud' , 'tmur' , 'LTG' , 'Y1G' , 'sub_Ret' , 'RMSE_PD' , 'exitflag'};%Variables to store

SimVar_Names_EIS_G = {'cg_yoy' , 'dg_yoy', 'rf' , 'EXRET' , 'exret' , 'ret_MS' , 'pd_ttm' , 'tmu' , 'tmud' , 'tmur' , 'LTG' , 'Y1G' , 'sub_ret'};%Variables to store

SimVar_Names_RE = {'cg_yoy' , 'dg_yoy', 'rf' , 'ret' , 'EXRET' , 'exret' , 'pd_ttm' , 'tmu' , 'tmud' , 'tmur'};%Variables to store

%% 2.2 EIS = 1

DATA_EIS_1 = Sim_EIS_1(l_q , l_q_burn , N_sim , nu_q , c_0 , node_n , SimVar_Names_EIS_1);

DATA_EIS_1_RE = Sim_EIS_1_RE(l_q , l_q_burn , N_sim , c_0 , SimVar_Names_RE);

%% 2.3 EIS ~= 1
DATA_EIS_G = Sim_EIS(l_q , l_q_burn , N_sim , nu_q , phi , psi , c_0 , SimVar_Names_EIS_G);

DATA_EIS_G_RE = Sim_EIS_RE(l_q , l_q_burn , N_sim , psi , c_0 , SimVar_Names_RE);


%% 2.4. Save simulated data for future use

% save Simulations_Large.mat SimVar_Names_EIS_1 DATA_EIS_1...
%                            SimVar_Names_EIS_G DATA_EIS_G...
%                            SimVar_Names_RE  DATA_EIS_1_RE DATA_EIS_G_RE ...
%      -v7.3 
               
%% ** 3. Table 5

%% 3.1 Calculate moments
Mom_EIS_1 = Mom_Annual(DATA_EIS_1 , type_avg , SimVar_Names_EIS_1 , MomVar_Names , Mom_RowNames , 'EIS 1');
Mom_EIS_1_RE = Mom_Annual(DATA_EIS_1_RE , type_avg , SimVar_Names_RE , MomVar_Names , Mom_RowNames , 'EIS 1 RE');


Mom_EIS_G = Mom_Annual(DATA_EIS_G , type_avg , SimVar_Names_EIS_G , MomVar_Names , Mom_RowNames , 'EIS General');
Mom_EIS_G_RE = Mom_Annual(DATA_EIS_G_RE , type_avg , SimVar_Names_RE , MomVar_Names , Mom_RowNames , 'EIS General RE');

%% 3.2 Output Table 5
Mom_list = {Mom_DATA , ...
            Mom_EIS_1 , Mom_EIS_1_RE ,...
            Mom_EIS_G , Mom_EIS_G_RE};

Mom_table = Mom_list{1};        
for i = 2 : length(Mom_list)
    [Mom_table , idx_temp] = outerjoin(Mom_table , Mom_list{i} , 'Keys','Row' ,'MergeKeys',1);
    [~ , reorder_temp] = sort(idx_temp);
    Mom_table = Mom_table(reorder_temp , :);
end


latextable2(table2array(Mom_table), 'Vert',Mom_RowNames,...
            'CustomRow',[0], 'CustomRowText',  ...
            {'& 1927-2019  & & \multicolumn{2}{c}{$\psi = 1$, $\phi = 1$} && \multicolumn{2}{c}{$\psi = 1.5$, $\phi = 0.99$}\\';...
            '\cline{4-5}\cline{7-8}\\' ;...
            '& & & Learning & RE & & Learning & RE'} ,...
            'format','%.2f','name', 'Results/Table_V.tex');
        
        
%% ** 4. Table 6

pred_horizon = [1 , 4 , 20];%Regression horizons to consider (in quarters)
mainRHS_6_name = {'tmud' , 'tmur' , 'tmu' ,'pd_ttm'};%Main RHS variables

%% 4.1 EIS = 1
% Extract simulated data
exret_EIS_1 = DATA_EIS_1{strcmp(SimVar_Names_EIS_1 , 'exret')};

Result_6_EIS_1 = cell(length(pred_horizon) , length(mainRHS_6_name));
for i = 1 : length(pred_horizon)%For a given horizon
    % moving sum
    exret_ms_EIS_1 = [nan(pred_horizon(i) - 1 , N_sim);... 
                      movsum(exret_EIS_1, pred_horizon(i) , 'Endpoints' , 'discard')];                  
    for j = 1 : length(mainRHS_6_name)%For each main RHS variable, e.g. expd
        mainRHS_full_temp = DATA_EIS_1{strcmp(SimVar_Names_EIS_1 , mainRHS_6_name{j})};%Extract main variables from the cell
        Result_sim_temp = nan(2 , N_sim);%Store results across simulations        
        for k = 1 : N_sim%For each simulation
            RHS_temp = mainRHS_full_temp(: , k);
            LHS_temp = exret_ms_EIS_1(: , k);
            
            [b_temp , ~ , ~ , Radj_temp] = olsgmm(LHS_temp(1 + pred_horizon(i) : end) ,...
                                                  [ones(l_q - pred_horizon(i) , 1) , RHS_temp(1 : end - pred_horizon(i))] ,...
                                                   0 , 1);
            Result_sim_temp(: , k) = [b_temp(2);Radj_temp];                                            
        end    
        
        Result_6_EIS_1{i , j} = mean(Result_sim_temp, 2);           
    end
end


%% 4.2 EIS ~= 1
% Extract simulated data
exret_EIS_G = DATA_EIS_G{strcmp(SimVar_Names_EIS_G , 'exret')};

Result_6_EIS_G = cell(length(pred_horizon) , length(mainRHS_6_name));
for i = 1 : length(pred_horizon)%For a given horizon
    % moving sum
    exret_ms_EIS_G = [nan(pred_horizon(i) - 1 , N_sim);... 
                      movsum(exret_EIS_G, pred_horizon(i) , 'Endpoints' , 'discard')];                  
    for j = 1 : length(mainRHS_6_name)%For each main RHS variable, e.g. expd
        mainRHS_full_temp = DATA_EIS_G{strcmp(SimVar_Names_EIS_G , mainRHS_6_name{j})};%Extract main variables from the cell
        Result_sim_temp = nan(2 , N_sim);%Store results across simulations        
        for k = 1 : N_sim%For each simulation
            RHS_temp = mainRHS_full_temp(: , k);
            LHS_temp = exret_ms_EIS_G(: , k);
            
            [b_temp , ~ , ~ , Radj_temp] = olsgmm(LHS_temp(1 + pred_horizon(i) : end) ,...
                                                  [ones(l_q - pred_horizon(i) , 1) , RHS_temp(1 : end - pred_horizon(i))] ,...
                                                   0 , 1);
            Result_sim_temp(: , k) = [b_temp(2);Radj_temp];                                            
        end    
        
        Result_6_EIS_G{i , j} = mean(Result_sim_temp, 2);           
    end
end

%% 4.3 Output Table 6

input_6.tableColLabels = {'(1)', '(2)' , '(3)' , '(4)' , '(5)' , '6'};

input_6.tableRowLabels = {'$\tilde{\mu}_d$' , '$R_{adj}^2$ ' , ... 
                          '$\tilde{\mu}_r$' , '$R_{adj}^2$ ' , ... 
                          '$p - d$'         , '$R_{adj}^2$ '};

input_6.data = [Result_6_EIS_1{1 , 1} , Result_6_EIS_1{2 , 1} , Result_6_EIS_1{3 , 1} , Result_6_EIS_G{1 , 1} , Result_6_EIS_G{2 , 1} , Result_6_EIS_G{3 , 1} ;...
                Result_6_EIS_1{1 , 2} , Result_6_EIS_1{2 , 2} , Result_6_EIS_1{3 , 2} , Result_6_EIS_G{1 , 2} , Result_6_EIS_G{2 , 2} , Result_6_EIS_G{3 , 2} ;...
                Result_6_EIS_1{1 , 4} , Result_6_EIS_1{2 , 4} , Result_6_EIS_1{3 , 4} , Result_6_EIS_G{1 , 4} , Result_6_EIS_G{2 , 4} , Result_6_EIS_G{3 , 4}];
                

input_6.dataFormatMode = 'row';
input_6.dataFormat = {'%.2f', 1 , '%.3f' , 1 , '%.2f', 1 , '%.3f' , 1 , '%.2f', 1 , '%.3f' , 1};   
                  
input_6.dataNanString = ' '; 
input_6.tableBorders = 0;

latex_Table_6 = latexTable(input_6 , 'Results/Table_VI.tex');

%% 5. Table 7

sub_Ret_EIS_1 = DATA_EIS_1{strcmp(SimVar_Names_EIS_1 , 'sub_Ret')};%Quarterly
Rf_EIS_1 = DATA_EIS_1{strcmp(SimVar_Names_EIS_1 , 'Rf')};%Quarterly
ret_MS_EIS_1 = DATA_EIS_1{strcmp(SimVar_Names_EIS_1 , 'ret_MS')};%Quarterly
Ret_MS_EIS_1 = exp(ret_MS_EIS_1) - 1;

mainRHS_7_name = {'tmud' , 'tmur'};%Main RHS variables
control_7 = ret_MS_EIS_1;

%% 5.1 Panel A
LHS_7A = (1 + sub_Ret_EIS_1).^4 - (1 + Rf_EIS_1).^4;

Result_7A = cell(length(mainRHS_7_name) + 1 , 1);
for j = 1 : length(mainRHS_7_name)
    mainRHS_full_temp = DATA_EIS_1{strcmp(SimVar_Names_EIS_1 , mainRHS_7_name{j})};
    
    Result_sim_temp = nan(3 * 2 , N_sim);
    for k = 1 : N_sim
        RHS_temp = [mainRHS_full_temp(: , k) , control_7(: , k)];
        for s = 1 : size(RHS_temp , 2)
            [b_temp , ~ , ~ , Radj_temp] = olsgmm(LHS_7A(: , k) , [ones(length(LHS_7A(: , k)) , 1) , RHS_temp(: , 1 : s)] , 0 , 1);
            Result_sim_temp((s - 1) * 3 + 1 : s * 3  , k) = [b_temp(2 : end);...
                                                             nan(size(RHS_temp , 2) - s , 1);...
                                                             Radj_temp];
        end
    end
    
    Result_7A{j} = reshape(mean(Result_sim_temp , 2) , 3 , []);
end

% Only lagged log returns
Result_7A_control = nan(3 , N_sim);
for k = 1 : N_sim
    [b_temp , ~ , ~ , Radj_temp] = olsgmm(LHS_7A(: , k) , [ones(length(LHS_7A(: , k)) , 1) , control_7(: , k)] , 0 , 1);
    Result_7A_control(: , k) = [NaN;
                                b_temp(2);...
                                Radj_temp];
end
    
Result_7A{end} = mean(Result_7A_control , 2);

%% 5.2 Panel B

LHS_7B = 1 + Ret_MS_EIS_1(5 : end , :) - (1 +sub_Ret_EIS_1(1 : end - 4 , :)).^4;

Result_7B = cell(length(mainRHS_7_name) + 1 , 1);
for j = 1 : length(mainRHS_7_name)
    mainRHS_full_temp = DATA_EIS_1{strcmp(SimVar_Names_EIS_1 , mainRHS_7_name{j})};
    
    Result_sim_temp = nan(3 * 2 , N_sim);
    for k = 1 : N_sim
        RHS_temp = [mainRHS_full_temp(1 : end - 4 , k) , control_7(1 : end - 4 , k)];
        for s = 1 : size(RHS_temp , 2)
            [b_temp , ~ , ~ , Radj_temp] = olsgmm(LHS_7B(: , k) , [ones(length(LHS_7B(: , k)) , 1) , RHS_temp(: , 1 : s)] , 0 , 1);
            Result_sim_temp((s - 1) * 3 + 1 : s * 3  , k) = [b_temp(2 : end);...
                                                             nan(size(RHS_temp , 2) - s , 1);...
                                                             Radj_temp];
        end
    end
    
    Result_7B{j} = reshape(mean(Result_sim_temp , 2) , 3 , []);
end

% Only lagged log returns
Result_7B_control = nan(3 , N_sim);
for k = 1 : N_sim
    [b_temp , ~ , ~ , Radj_temp] = olsgmm(LHS_7B(: , k) , [ones(length(LHS_7B(: , k)) , 1) , control_7(1 : end - 4 , k)] , 0 , 1);
    Result_7B_control(: , k) = [NaN;
                                b_temp(2);...
                                Radj_temp];
end
    
Result_7B{end} = mean(Result_7B_control , 2);

%% 5.3 Output Table 7
input_7.tableColLabels = {'(1)', '(2)' , '(3)' , '(4)' , '(5)'};

input_7.tableRowLabels = repmat({'$\tilde{\mu}$' , ... 
                                 '$r_{t-3,t}$' , ... 
                                 'Adj. $R^2$ '} , 1 , 2);

input_7.data = [Result_7A{1}(: , 1) , Result_7A{end}(: , 1) , Result_7A{1}(: , 2) , Result_7A{2};...
                Result_7B{1}(: , 1) , Result_7B{end}(: , 1) , Result_7B{1}(: , 2) , Result_7B{2}];
                

input_7.dataFormatMode = 'row';
input_7.dataFormat = {'%.3f', 2 , '%.3f' , 1 , '%.3f', 2 , '%.3f' , 1};   
                  
input_7.dataNanString = ' '; 
input_7.tableBorders = 0;

latex_Table_7 = latexTable(input_7 , 'Results/Table_VII.tex');

%% 6. Footnote 22

% Conditional variance on the tmu_d as in Table VIII
N_cond_sim = 1000;
N_sim_2 = 100;

SimVar_Names_EIS_1_CV = {'cg_yoy' , 'dg_yoy', 'rf' , 'Rf' , 'Ret' , 'ret' , 'ret_MS' , 'EXRET' , 'exret' , 'pd_ttm' , 'tmu' , 'tmud' , 'tmur' , 'condVar'};%Variables to store

DATA_EIS_1_CV = Sim_EIS_1_CV(l_q , l_q_burn , N_sim_2 , N_cond_sim , nu_q , c_0 , SimVar_Names_EIS_1_CV);

LHS_CV = DATA_EIS_1_CV{strcmp(SimVar_Names_EIS_1_CV , 'condVar')};
RHS_CV_tmud = DATA_EIS_1_CV{strcmp(SimVar_Names_EIS_1_CV , 'tmud')};
RHS_CV_pd = DATA_EIS_1_CV{strcmp(SimVar_Names_EIS_1_CV , 'pd_ttm')}/100;
RHS_CV_rf = DATA_EIS_1_CV{strcmp(SimVar_Names_EIS_1_CV , 'rf')};

result_condVar_tmud  = nan(9 , N_sim_2); % Only use tmud
result_condVar_pd    = nan(9 , N_sim_2); % Only use pd
result_condVar_pd_rf = nan(9 , N_sim_2); % Use pd and rf

for i = 1 : N_sim_2
    [b_temp , ~ , p_temp , Radj_temp] = olsnan_gmm_EWC(LHS_CV(: , i) , RHS_CV_tmud(: , i) , 1 , 1);
    result_condVar_tmud([7 1], i) = b_temp;  
    result_condVar_tmud([8 2], i) = p_temp;  
    result_condVar_tmud(9, i) = Radj_temp;  
    
    [b_temp , ~ , p_temp , Radj_temp] = olsnan_gmm_EWC(LHS_CV(: , i) , RHS_CV_pd(: , i) , 1 , 1);
    result_condVar_pd([7 3], i) = b_temp;  
    result_condVar_pd([8 4], i) = p_temp;  
    result_condVar_pd(9, i) = Radj_temp;     
    
    [b_temp , ~ , p_temp , Radj_temp] = olsnan_EWC(LHS_CV(: , i) , [RHS_CV_pd(: , i) , RHS_CV_rf(: , i)] , 1 , 1);
    result_condVar_pd_rf([7 3 5], i) = b_temp;  
    result_condVar_pd_rf([8 4 6], i) = p_temp;  
    result_condVar_pd_rf(9, i) = Radj_temp;   
end

disp('-------------------------------')
disp('Mean of conditional variance')
disp(mean(mean(LHS_CV)) * 100)
disp('Median s.d. of conditional variance')
disp(mean(std(LHS_CV)) * 100)
disp('Mean coefficients of conditional variance on tmud')
disp(mean(result_condVar_tmud(1 , :)))
disp('Change of conditional variance with one s.d. of tmud')
disp(mean(result_condVar_tmud(1 , :)) * 1.31/4)

%% 7. Figure 3

%% 7.1 Get statistics
OOS_Var = {'pd_ttm' , 'tmud'};
OOS_Var_sign = [-1 , -1];
OOS_horizon = 4;
OOS_burn = 80;

[IS_path , OOS_path] = OOS_EIS_1(DATA_EIS_1 , SimVar_Names_EIS_1 , OOS_Var , OOS_Var_sign ,OOS_horizon , OOS_burn);

%% 7.2 Output Figure 3

figure
% pd
subplot(2 , 1 , 1)
plot((OOS_horizon + 1 : l_q)' , IS_path{1} , 'b' , (OOS_horizon + OOS_burn : l_q)' , OOS_path{1 , 1} , 'r' , (OOS_horizon + OOS_burn : l_q)' , OOS_path{1 , 2} , '-.k' , 'LineWidth' , 2)
l1 = legend('IS' , 'OOS (Unconstrained)' , 'OOS (Constrained)','Location','northwest','AutoUpdate','off');
l1.FontSize = 10;
title('$p - d$','interpreter','latex','FontSize',12)
xL = xlim;
ylim([-0.2 0.5])
line(xL, [0 0],'Color',[0, 0, 0, 0.2]);  
xlabel('Period','FontSize',10)
ylabel('Cumulative SPE Difference','FontSize',10) 
% tmud
subplot(2 , 1 , 2)
plot((OOS_horizon + 1 : l_q)' , IS_path{2} , 'b' , (OOS_horizon + OOS_burn : l_q)' , OOS_path{2 , 1} , 'r' , (OOS_horizon + OOS_burn : l_q)' , OOS_path{2 , 2} , '-.k' , 'LineWidth' , 2)
l1 = legend('IS' , 'OOS (Unconstrained)' , 'OOS (Constrained)','Location','northwest','AutoUpdate','off');
l1.FontSize = 10;
title('$\tilde{\mu}_d$','interpreter','latex','FontSize',12)
xL = xlim;
ylim([-0.2 0.5])
line(xL, [0 0],'Color',[0, 0, 0, 0.2]);  
xlabel('Period','FontSize',10)
ylabel('Cumulative SPE Difference','FontSize',10) 
set(gca,'TickLabelInterpreter', 'latex')
saveas(gcf,'Figures/Figure3','epsc')
savetightfigure('Figures/Figure3')

%% 8. Results in Section 3.4 (Long-term Subjective Dividend Growth Expectations)

%% 8.1 EIS = 1
LTG_EIS1_sim = nan(12 , N_sim); % Corresponds to
                             % 1) tmud
                             % 2) tmud and one-year lagged return
                             % 3) tmur
                             % 4) tmur and one-year lagged return
                             % 5) Only lagged one-year return
                             
LTG_EIS1 = DATA_EIS_1{strcmp(SimVar_Names_EIS_1 , 'LTG')};
RetMS_EIS1 = exp(DATA_EIS_1{strcmp(SimVar_Names_EIS_1 , 'ret_MS')}) - 1;
tmud_EIS1 = DATA_EIS_1{strcmp(SimVar_Names_EIS_1 , 'tmud')};
tmur_EIS1 = DATA_EIS_1{strcmp(SimVar_Names_EIS_1 , 'tmur')};
for i = 1 : N_sim
    % Only tmud
    [b_1_temp , ~ , ~ , Radj_1_temp] = olsgmm(LTG_EIS1(: , i) , [ones(l_q , 1) , tmud_EIS1(: , i)] , 0 , 1); 
    LTG_EIS1_sim(1 , i) = b_1_temp(2);
    LTG_EIS1_sim(2 , i) = Radj_1_temp;
    % tmud and one-year lagged return
    [b_2_temp , ~ , ~ , Radj_2_temp] = olsgmm(LTG_EIS1(: , i) , [ones(l_q , 1) , tmud_EIS1(: , i) , RetMS_EIS1(: , i)] , 0 , 1);
    LTG_EIS1_sim(3 : 4 , i) = b_2_temp(2 : 3);
    LTG_EIS1_sim(5 , i) = Radj_2_temp;

    % Only tmur
    [b_3_temp , ~ , ~ , Radj_3_temp] = olsgmm(LTG_EIS1(: , i) , [ones(l_q , 1) , tmur_EIS1(: , i)] , 0 , 1);
    LTG_EIS1_sim(6 , i) = b_3_temp(2);
    LTG_EIS1_sim(7 , i) = Radj_3_temp;
    
    % tmur and one-year lagged return
    [b_4_temp , ~ , ~ , Radj_4_temp] = olsgmm(LTG_EIS1(: , i) , [ones(l_q , 1) , tmur_EIS1(: , i) , RetMS_EIS1(: , i)] , 0 , 1);
    LTG_EIS1_sim(8 : 9 , i) = b_4_temp(2 : 3);
    LTG_EIS1_sim(10 , i) = Radj_4_temp;
    
    % one-year lagged return
    [b_5_temp , ~ , ~ , Radj_5_temp] = olsgmm(LTG_EIS1(: , i) , [ones(l_q , 1) , RetMS_EIS1(: , i)] , 0 , 1);
    LTG_EIS1_sim(11 , i) = b_5_temp(2);
    LTG_EIS1_sim(12 , i) = Radj_5_temp;
end

LTG_EIS1_result = mean(LTG_EIS1_sim , 2);

%% 8.2 EIS ~= 1
LTG_EISG_sim = nan(12 , N_sim); % Corresponds to
                             % 1) tmud
                             % 2) tmud and one-year lagged return
                             % 3) tmur
                             % 4) tmur and one-year lagged return
                             % 5) Only lagged one-year return
                             
LTG_EISG = DATA_EIS_G{strcmp(SimVar_Names_EIS_G , 'LTG')};
RetMS_EISG = exp(DATA_EIS_G{strcmp(SimVar_Names_EIS_G , 'ret_MS')}) - 1;
tmud_EISG = DATA_EIS_G{strcmp(SimVar_Names_EIS_G , 'tmud')};
tmur_EISG = DATA_EIS_G{strcmp(SimVar_Names_EIS_G , 'tmur')};
for i = 1 : N_sim
    % Only tmud
    [b_1_temp , ~ , ~ , Radj_1_temp] = olsgmm(LTG_EISG(: , i) , [ones(l_q , 1) , tmud_EISG(: , i)] , 0 , 1); 
    LTG_EISG_sim(1 , i) = b_1_temp(2);
    LTG_EISG_sim(2 , i) = Radj_1_temp;
    % tmud and one-year lagged return
    [b_2_temp , ~ , ~ , Radj_2_temp] = olsgmm(LTG_EISG(: , i) , [ones(l_q , 1) , tmud_EISG(: , i) , RetMS_EISG(: , i)] , 0 , 1);
    LTG_EISG_sim(3 : 4 , i) = b_2_temp(2 : 3);
    LTG_EISG_sim(5 , i) = Radj_2_temp;

    % Only tmur
    [b_3_temp , ~ , ~ , Radj_3_temp] = olsgmm(LTG_EISG(: , i) , [ones(l_q , 1) , tmur_EISG(: , i)] , 0 , 1);
    LTG_EISG_sim(6 , i) = b_3_temp(2);
    LTG_EISG_sim(7 , i) = Radj_3_temp;
    
    % tmur and one-year lagged return
    [b_4_temp , ~ , ~ , Radj_4_temp] = olsgmm(LTG_EISG(: , i) , [ones(l_q , 1) , tmur_EISG(: , i) , RetMS_EISG(: , i)] , 0 , 1);
    LTG_EISG_sim(8 : 9 , i) = b_4_temp(2 : 3);
    LTG_EISG_sim(10 , i) = Radj_4_temp;
    
    % one-year lagged return
    [b_5_temp , ~ , ~ , Radj_5_temp] = olsgmm(LTG_EISG(: , i) , [ones(l_q , 1) , RetMS_EISG(: , i)] , 0 , 1);
    LTG_EISG_sim(11 , i) = b_5_temp(2);
    LTG_EISG_sim(12 , i) = Radj_5_temp;
end

LTG_EISG_result = mean(LTG_EISG_sim , 2);

%% 8.3 Output Table LTG

input_LTG.tableColLabels = {'(1)', '(2)' , '(3)' , '(4)' , '(5)'};

input_LTG.tableRowLabels = repmat({'$\tilde{\mu}_d$' , ... 
                                 '$\tilde{\mu}_r$' , ... 
                                 '$r_{t-3,t}$' , ... 
                                 'Adj. $R^2$ '} , 1 , 2);

input_LTG.data = [LTG_EIS1_result(1) , NaN                , LTG_EIS1_result(3) , NaN                , NaN;...
                NaN                , NaN                  , NaN                , LTG_EIS1_result(6) , LTG_EIS1_result(8);...
                NaN                , LTG_EIS1_result(11)  , LTG_EIS1_result(4) , NaN                , LTG_EIS1_result(9);...
                LTG_EIS1_result(2) , LTG_EIS1_result(12)  , LTG_EIS1_result(5) , LTG_EIS1_result(7) , LTG_EIS1_result(10);...
                LTG_EISG_result(1) , NaN                  , LTG_EISG_result(3) , NaN                , NaN;...
                NaN                , NaN                  , NaN                , LTG_EISG_result(6) , LTG_EISG_result(8);...
                NaN                , LTG_EISG_result(11)  , LTG_EISG_result(4) , NaN                , LTG_EISG_result(9);...
                LTG_EISG_result(2) , LTG_EISG_result(12)  , LTG_EISG_result(5) , LTG_EISG_result(7) , LTG_EISG_result(10)];
                

input_LTG.dataFormatMode = 'row';
input_LTG.dataFormat = {'%.2f', 3 , '%.3f' , 1 , '%.2f', 3 , '%.3f' , 1};   
                  
input_LTG.dataNanString = ' '; 
input_LTG.tableBorders = 0;

latex_Table_LTG = latexTable(input_LTG , 'Results/Table_LTG.tex');

%% 9. Table D.1

disp('fmincon exitflag')
disp(DATA_EIS_1{strcmp(SimVar_Names_EIS_1 , 'exitflag')})

RMSE_PD_EIS_1 = DATA_EIS_1{strcmp(SimVar_Names_EIS_1 , 'RMSE_PD')};
disp('------------------------')
disp('RMSE of P/D ratio from Projection')
disp('Mean')
disp(mean(RMSE_PD_EIS_1))
disp('Median')
disp(median(RMSE_PD_EIS_1))
disp('Max')
disp(max(RMSE_PD_EIS_1))

I_param = param_Data;
tmu_EIS_1 = DATA_EIS_1{strcmp(SimVar_Names_EIS_1 , 'tmu')};
disp('------------------------')
disp('Grid for mu')
disp([I_param.mu - 4 * mean(std(tmu_EIS_1),2) , I_param.mu + 4 * mean(std(tmu_EIS_1),2)])
disp('Grid for d - c')
disp([I_param.d_c_mean - 4 * I_param.d_c_mean_sig , I_param.d_c_mean + 4 * I_param.d_c_mean_sig])

